---

Drawing maps

---

Here we will draw maps given a set of coordinates, using free services, and in several different styles. We will also use geolocation to find coordinates given the name of a place.

Import libraries

Install them if necessary:

#! conda install -y pandas matplotlib plotly geopy cartopy

Now import:

import plotly.express as px
import pandas as pd
import geopy
import cartopy
import matplotlib.pyplot as plt

df = pd.DataFrame([
  {'place':"AMNH"                 , 'lat':40.7813295, 'lon':-73.9741718, 'height_m':48,  'type':"Building"}, # made up height!
  {'place':"Statue of Liberty"    , 'lat':40.6893490, 'lon':-74.0444827, 'height_m':93,  'type':"Monument"},
  {'place':"Empire State Building", 'lat':40.7484445, 'lon':-73.9878584, 'height_m':443, 'type':"Building"},
  {'place':"Wasghington Statue"   , 'lat':40.7073993, 'lon':-74.0123975, 'height_m':18,  'type':"Monument"}, # also made up
])
df

	place	lat	lon	height_m	type
0	AMNH	40.781329	-73.974172	48	Building
1	Statue of Liberty	40.689349	-74.044483	93	Monument
2	Empire State Building	40.748444	-73.987858	443	Building
3	Wasghington Statue	40.707399	-74.012398	18	Monument

Plot a zoomable USGS-based map with Plotly

fig = px.scatter_mapbox(df,
                        lat='lat',
                        lon='lon',
                        color='type',
                        hover_name='place',
                        text='place',
                        size='height_m',
                        zoom=10,
                        color_discrete_sequence=px.colors.qualitative.Alphabet_r
                       )

fig.layout.title          = "NYC Landmarks"
fig.layout.mapbox.style   = 'white-bg'
fig.layout.mapbox.layers  = [
    {
        'below'             : 'traces',
        'sourcetype'        : 'raster',
        'sourceattribution' : "",#"United States Geological Survey",
        'source'            : ['https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/tile/{z}/{y}/{x}']
    }
]

fig.show()

Open Street Maps

fig = px.scatter_mapbox(df, lat='lat', lon='lon', hover_name='place', hover_data=['place'],
                        zoom=10, height=300, size='height_m',
                        color_discrete_sequence=['green', 'firebrick'], color='type')
fig.update_layout(mapbox_style="open-street-map")
fig.update_layout(margin={'r':0,'t':0,'l':0,'b':0})
fig.show()

Using the stamen-terrain variation

fig = px.scatter_mapbox(df, lat='lat', lon='lon', hover_name='place', hover_data=['place'],
                     color='type', zoom=10, height=500, size='height_m')

fig.update_layout(mapbox_style="stamen-terrain",# mapbox_zoom=4, mapbox_center_lat = 41,
    margin={"r":0,"t":0,"l":0,"b":0})

fig.show()

Export the plotly map to an HTML file

# Constants to hide the Plotly navigation bar at the top, or just the logo
PLOTLY_CONFIG_NOBAR  = {'displayModeBar' : False}
PLOTLY_CONFIG_NOLOGO = {'displaylogo'    : False}
# Export to HTML
fig.write_html('my_maps.html', config=PLOTLY_CONFIG_NOLOGO)

Note: exporting a Plotly map to PNG is a lot more difficult, and I do not cover it here (yet).

Static maps using Matplotlib & Cartopy

Add layers to a basic map

Different layers can be added with add_feature. Below, we add LAND and OCEAN.

plt.figure(figsize=(18,8))
axes = plt.axes(projection=cartopy.crs.PlateCarree())
axes.add_feature(cartopy.feature.COASTLINE, linewidth=0.4)
#axes.add_feature(cartopy.feature.BORDERS, linewidth=0.6)
axes.add_feature(cartopy.feature.LAND)
axes.add_feature(cartopy.feature.OCEAN)

for i in df.index:
  axes.plot(df.iloc[i]['lon'], df.iloc[i]['lat'], 'b.')
  
plt.axis(xmin=-75, xmax=-73, ymin=40.5, ymax=41)

(-75.0, -73.0, 40.5, 41.0)

---

Not super pretty, but it works.

A natural-surface style

ax = plt.axes(projection=cartopy.crs.Mollweide())
ax.stock_img()
plt.show()

This looks very nice for continental distances, but if we try to use it for a dataset of nearby points like the one for NYC that we made above, it pixelates so horrendously it becomes useless:

plt.figure(figsize=(18,8))
axes = plt.axes(projection=cartopy.crs.PlateCarree())
axes.stock_img()

for i in df.index:
  axes.plot(df.iloc[i]['lon'], df.iloc[i]['lat'], 'b.')
  
plt.axis(xmin=-75, xmax=-73, ymin=40.5, ymax=41)

(-75.0, -73.0, 40.5, 41.0)

With this kind of distance, this map style is not useful. We could try to zoom out by changing the x and y limits, but then we lose the points we were trying to observe:

plt.figure(figsize=(18,8))
axes = plt.axes(projection=cartopy.crs.PlateCarree())
axes.stock_img()

for i in df.index:
  axes.plot(df.iloc[i]['lon'], df.iloc[i]['lat'], 'b.')
  
plt.axis(xmin=-85, xmax=-63, ymin=30.5, ymax=51)

(-85.0, -63.0, 30.5, 51.0)

P.S.: Pandas provides much more efficient ways of plotting the points than the for loop used above, but this is not a tutorial about Pandas.

Finding locations by name with a geocoder

Just like map providers, there are also many geolocators. You're no doubt familiar with Google Maps. We could access their geolocation system through Python but, unlike the on the Google Maps website, the programmatic one isn't free. There are options, however.

We imported geopy above, which has a geocoder called Nominatim that can find locations by name.

Let's try it with a few trivial and less-than-trivial names:

• Caracas, Antananarivo and Jakarta should be easy, as there's only one major one of each those.
• Nueva York is not written in English. Can the system find it?
• Minga is the name given to Munich only in the local Bavarian dialect, and even most Germans don't know that; let's see how the searcher fares.
• There are obvious choices for Rome, Wellington and Manchester in Italy, New Zealand and England, but there are several identically named cities across North America; let's see if the geolocator gets confused.
• Springfield is a bit of a mean one, since there is at least one of those in every state of the USA, plus several more in the UK, Canada and Australia. Which one will the system pick?
• The Big Easy is a nickname often given to New Orleans. Let's see if the system picks it up.
• Taj Mahal, Maracana and The Great Sphinx are not cities, but they are famous landmarks that should, at least in principle, be easy to find.

geolocator = geopy.geocoders.Nominatim(user_agent='learning-maps')

cities = ['Caracas', 'Jakarta', 'Nueva York', 'Minga', 'Rome', 'Manchester',
          'Springfield', 'The Big Easy', 'Wellington', 'Antananarivo',
          'Taj Mahal', 'Maracana', 'The Great Sphinx']

# An empty directory, to be filled below with the locator for each city
locations = {}

for city in cities:
  print(f"Geolocating '{city}'...")
  locations[city] = geolocator.geocode(city)
  print("Found at:", locations[city], "\n")

Geolocating 'Caracas'...
Found at: Caracas, Parroquia Catedral, Municipio Libertador, Distrito Capital, Venezuela 

Geolocating 'Jakarta'...
Found at: Daerah Khusus Ibukota Jakarta, Indonesia 

Geolocating 'Nueva York'...
Found at: New York, United States 

Geolocating 'Minga'...
Found at: München, Bayern, Deutschland 

Geolocating 'Rome'...
Found at: Roma, Roma Capitale, Lazio, Italia 

Geolocating 'Manchester'...
Found at: Manchester, Greater Manchester, North West England, England, United Kingdom 

Geolocating 'Springfield'...
Found at: Springfield, Sangamon County, Illinois, United States 

Geolocating 'The Big Easy'...
Found at: The Big Easy, Dalbeattie, Dumfries and Galloway, Scotland, DG5 4QR, United Kingdom 

Geolocating 'Wellington'...
Found at: Wellington, Wellington City, Wellington, 6011, New Zealand / Aotearoa 

Geolocating 'Antananarivo'...
Found at: Antananarivo, Analamanga, Province d’Antananarivo, 101, Madagasikara 

Geolocating 'Taj Mahal'...
Found at: Taj Mahal, Taj Mahal Internal Path, Taj Ganj, Agra, Uttar Pradesh, 282001, India 

Geolocating 'Maracana'...
Found at: Maracanã, Avenida Presidente Castelo Branco, Maracanã, Zona Norte do Rio de Janeiro, Rio de Janeiro, Região Geográfica Imediata do Rio de Janeiro, Região Metropolitana do Rio de Janeiro, Região Geográfica Intermediária do Rio de Janeiro, Rio de Janeiro, Região Sudeste, 20271-130, Brasil 

Geolocating 'The Great Sphinx'...
Found at: Great Sphinx, Hill Country Retreat, Bexar County, Texas, 78253, United States 

---

Mostly great in my case, except for The Big Easy and The Great Sphinx, which were complete failures (it may be different when you run it). We should probably go back and fix those (e.g., add "of Giza" to the sphinx). Note also that the Springfield in Illinois was favoured by the geolocator, which perhaps makes sense as that one is a State capital.

In fact, we have the coordinates:

locations['Nueva York'].latitude, locations['Nueva York'].longitude

(40.7127281, -74.0060152)

This would be a much better use case for our nice National-Geographic-looking map:

plt.figure(figsize=(18,8))
axes = plt.axes(projection=cartopy.crs.PlateCarree())
axes.stock_img()

for city, location in locations.items():
  # Put a point on the map for each city
  axes.plot(location.longitude, location.latitude, 'b.') # 'b.' means "blue dots"
  # Add the name at the same position
  axes.text(location.longitude, location.latitude, city)

Well, the positioning and sizing of the texts will need some work, but this is not a mabplotlib tutorial, so this will do for now. Also remember that we need to fix the Sphinx and New Orleans, both of which are utterly wrong.